# THIS SCRIPT CREATES TABLE 3 and TABLE A.4

#clear environment 
rm(list = ls())

library(readr)
library(tidyr)
library(haven)

require(tigris) #
require(tidycensus) #
require(sp)
require(rgdal) #
require(sf) #
require(raster)
require(stargazer) #
require(readxl)
require(haven)
require(zoo) #== for na.locf
require(snow) #
require(doSNOW) #
require(parallel)
require(foreach)
require(iterators)
require(foreign)
require(abind)
library(dplyr)
library(Rfast)

#set directories
WORK.DIRECTORY = "Z:/user/rih4/Solar Siting/JAERE Code Repository"
WORK.OUT = file.path(WORK.DIRECTORY,"Data JAERE/intermediate_output")
PLOT.OUT = file.path(WORK.DIRECTORY, "Data JAERE/final_output")
dir.create(WORK.OUT);dir.create(PLOT.OUT, recursive=T)

reallocation.interstate <- read_csv(file.path(WORK.OUT,"Simulation_Results_Country.csv")) %>% 
  mutate(zip = sprintf("%05d", zipcode)) %>%
  dplyr::select(-c(zipcode)) %>%
  relocate(zip) %>%
  replace_na(list(num_systems=0)) %>%
  mutate(total_damages_avoided_original = total_damages*num_systems,
         total_damages_avoided_interstate_env = total_damages*new_num_systems,
         total_damages_avoided_interstate_env30 = total_damages*new_num_systems30)

reallocation.intrastate <- read_csv(file.path(WORK.OUT,"Simulation_Results_Intrastate.csv")) %>% 
  mutate(zip = sprintf("%05d", zipcode)) %>%
  dplyr::select(-c(zipcode)) %>%
  relocate(zip) %>%
  replace_na(list(num_systems=0)) %>%
  mutate(total_damages_avoided_original = total_damages*num_systems,
         total_damages_avoided_intrastate_env = total_damages*state_new_num_systems,
         total_damages_avoided_intrastate_env30 = total_damages*state_new_num_systems30) %>%
  select(-c(NAME,state,POPULATION,POP_SQMI,pcac,INT,total_damages,num_systems,total_damages_avoided_original))

total.damages.env.all <- left_join(reallocation.interstate,reallocation.intrastate,by="zip") %>%
  summarize(original = sum(total_damages_avoided_original),
            interstate = sum(total_damages_avoided_interstate_env),
            interstate30 = sum(total_damages_avoided_interstate_env30),
            intrastate = sum(total_damages_avoided_intrastate_env),
            intrastate30 = sum(total_damages_avoided_intrastate_env30)) %>%
  mutate(across(everything(),~round(.x/1000000,digits=1)))

number.zips.env.all <- left_join(reallocation.interstate,reallocation.intrastate,by="zip") %>%
  summarize(original = sum(num_systems>0),
            interstate = sum(new_num_systems>0),
            interstate30 = sum(new_num_systems30>0),
            intrastate = sum(state_new_num_systems>0),
            intrastate30 = sum(state_new_num_systems30>0))

mean.damages.env.all <- left_join(reallocation.interstate,reallocation.intrastate,by="zip") %>%
  summarize(original = mean(total_damages[num_systems>0]),
            interstate = mean(total_damages[new_num_systems>0]),
            interstate30 = mean(total_damages[new_num_systems30>0]),
            intrastate = mean(total_damages[state_new_num_systems>0]),
            intrastate30 = mean(total_damages[state_new_num_systems30>0])) %>%
  mutate(across(everything(),round,digits=1))

table3 <- total.damages.env.all %>%
  bind_rows(number.zips.env.all,mean.damages.env.all) %>%
  mutate(Var = c("total avoided damages","number of zips","mean avoided damages")) %>%
  relocate(Var)
write_csv(table3,file.path(PLOT.OUT,"table3.csv"))

reallocation.interstate.env_energy <- read_csv(file.path(WORK.OUT,"Simulation_Results_Country_energy_plus_env.csv")) %>% 
  mutate(zip = sprintf("%05d", zipcode)) %>%
  dplyr::select(-c(zipcode)) %>%
  relocate(zip) %>%
  replace_na(list(num_systems=0)) %>%
  mutate(total_damages_avoided_original = total_damages*num_systems,
         total_damages_avoided_interstate_env_energy = total_damages*new_num_systems,
         total_damages_avoided_interstate_env_energy30 = total_damages*new_num_systems30) %>%
  rename(total_damages_plus_energy = total_damages)

reallocation.intrastate.env_energy <- read_csv(file.path(WORK.OUT,"Simulation_Results_Intrastate_energy_plus_env.csv")) %>% 
  mutate(zip = sprintf("%05d", zipcode)) %>%
  dplyr::select(-c(zipcode)) %>%
  relocate(zip) %>%
  replace_na(list(num_systems=0)) %>%
  mutate(total_damages_avoided_original = total_damages*num_systems,
         total_damages_avoided_intrastate_env_energy = total_damages*state_new_num_systems,
         total_damages_avoided_intrastate_env_energy30 = total_damages*state_new_num_systems30) %>%
  select(-c(NAME,state,POPULATION,POP_SQMI,pcac,INT,total_damages,num_systems,total_damages_avoided_original))

total.damages.env.energy.all <- left_join(reallocation.interstate.env_energy,reallocation.intrastate.env_energy,by="zip") %>%
  summarize(original = sum(total_damages_avoided_original),
            interstate = sum(total_damages_avoided_interstate_env_energy),
            interstate30 = sum(total_damages_avoided_interstate_env_energy30),
            intrastate = sum(total_damages_avoided_intrastate_env_energy),
            intrastate30 = sum(total_damages_avoided_intrastate_env_energy30)) %>%
  mutate(across(everything(),~round(.x/1000000,digits=1)))

number.zips.env.energy.all <- left_join(reallocation.interstate.env_energy,reallocation.intrastate.env_energy,by="zip") %>%
  summarize(original = sum(num_systems>0),
            interstate = sum(new_num_systems>0),
            interstate30 = sum(new_num_systems30>0),
            intrastate = sum(state_new_num_systems>0),
            intrastate30 = sum(state_new_num_systems30>0))

mean.damages.env.energy.all <- left_join(reallocation.interstate.env_energy,reallocation.intrastate.env_energy,by="zip") %>%
  summarize(original = mean(total_damages_plus_energy[num_systems>0]),
            interstate = mean(total_damages_plus_energy[new_num_systems>0]),
            interstate30 = mean(total_damages_plus_energy[new_num_systems30>0]),
            intrastate = mean(total_damages_plus_energy[state_new_num_systems>0]),
            intrastate30 = mean(total_damages_plus_energy[state_new_num_systems30>0])) %>%
  mutate(across(everything(),round,digits=1))

tableA4 <- total.damages.env.energy.all %>%
  bind_rows(number.zips.env.energy.all,mean.damages.env.energy.all) %>%
  mutate(Var = c("total avoided damages plus energy value","number of zips","mean avoided damages")) %>%
  relocate(Var)
write_csv(tableA4,file.path(PLOT.OUT,"tableA4.csv"))
